home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / hobby / ast44src.zip / MATRIX.C < prev    next >
Encoding:
C/C++ Source or Header  |  1995-02-11  |  21.7 KB  |  731 lines
  1. /*
  2. ** Astrolog (Version 4.40) File: matrix.c
  3. **
  4. ** IMPORTANT NOTICE: The graphics database and chart display routines
  5. ** used in this program are Copyright (C) 1991-1995 by Walter D. Pullen
  6. ** (astara@u.washington.edu). Permission is granted to freely use and
  7. ** distribute these routines provided one doesn't sell, restrict, or
  8. ** profit from them in any way. Modification is allowed provided these
  9. ** notices remain with any altered or edited versions of the program.
  10. **
  11. ** The main planetary calculation routines used in this program have
  12. ** been Copyrighted and the core of this program is basically a
  13. ** conversion to C of the routines created by James Neely as listed in
  14. ** Michael Erlewine's 'Manual of Computer Programming for Astrologers',
  15. ** available from Matrix Software. The copyright gives us permission to
  16. ** use the routines for personal use but not to sell them or profit from
  17. ** them in any way.
  18. **
  19. ** The PostScript code within the core graphics routines are programmed
  20. ** and Copyright (C) 1992-1993 by Brian D. Willoughby
  21. ** (brianw@sounds.wa.com). Conditions are identical to those above.
  22. **
  23. ** The extended accurate ephemeris databases and formulas are from the
  24. ** calculation routines in the program "Placalc" and are programmed and
  25. ** Copyright (C) 1989,1991,1993 by Astrodienst AG and Alois Treindl
  26. ** (alois@azur.ch). The use of that source code is subject to
  27. ** regulations made by Astrodienst Zurich, and the code is not in the
  28. ** public domain. This copyright notice must not be changed or removed
  29. ** by any user of this program.
  30. **
  31. ** Initial programming 8/28,30, 9/10,13,16,20,23, 10/3,6,7, 11/7,10,21/1991.
  32. ** X Window graphics initially programmed 10/23-29/1991.
  33. ** PostScript graphics initially programmed 11/29-30/1992.
  34. ** Last code change made 1/29/1995.
  35. */
  36.  
  37. #include "astrolog.h"
  38.  
  39.  
  40. #ifdef MATRIX
  41. real MC, Asc, RA, OB;
  42.  
  43.  
  44. /*
  45. ******************************************************************************
  46. ** Assorted Calculations.
  47. ******************************************************************************
  48. */
  49.  
  50. /* Given a month, day, and year, convert it into a single Julian day value, */
  51. /* i.e. the number of days passed since a fixed reference date.             */
  52.  
  53. long MdyToJulian(mon, day, yea)
  54. int mon, day, yea;
  55. {
  56. #ifndef PLACALC
  57.   long im, j;
  58.  
  59.   im = 12*((long)yea+4800)+(long)mon-3;
  60.   j = (2*(im%12) + 7 + 365*im)/12;
  61.   j += (long)day + im/48 - 32083;
  62.   if (j > 2299171)                   /* Take care of dates in */
  63.     j += im/4800 - im/1200 + 38;     /* Gregorian calendar.   */
  64.   return j;
  65. #else
  66.   int fGreg = fTrue;
  67.  
  68.   if (yea < yeaJ2G || (yea == yeaJ2G &&
  69.     (mon < monJ2G || (mon == monJ2G && day < 15))))
  70.     fGreg = fFalse;
  71.   return (long)RFloor(julday(mon, day, yea, 12.0, fGreg)+rRound);
  72. #endif
  73. }
  74.  
  75.  
  76. /* Like above but return a fractional Julian time given the extra info. */
  77.  
  78. real MdytszToJulian(mon, day, yea, tim, dst, zon)
  79. int mon, day, yea;
  80. real tim, dst, zon;
  81. {
  82.   return (real)MdyToJulian(mon, day, yea) +
  83.     (DecToDeg(tim) + DecToDeg(zon) - DecToDeg(dst)) / 24.0;
  84. }
  85.  
  86.  
  87. /* Take a Julian day value, and convert it back into the corresponding */
  88. /* month, day, and year.                                               */
  89.  
  90. void JulianToMdy(JD, mon, day, yea)
  91. real JD;
  92. int *mon, *day, *yea;
  93. {
  94. #ifndef PLACALC
  95.   long L, N, IT, JT, K, IK;
  96.  
  97.   L  = (long)RFloor(JD+rRound)+68569L;
  98.   N  = Dvd(4L*L, 146097L);
  99.   L  -= Dvd(146097L*N + 3L, 4L);
  100.   IT = Dvd(4000L*(L+1L), 1461001L);
  101.   L  -= Dvd(1461L*IT, 4L) - 31L;
  102.   JT = Dvd(80L*L, 2447L);
  103.   K  = L-Dvd(2447L*JT, 80L);
  104.   L  = Dvd(JT, 11L);
  105.   JT += 2L - 12L*L;
  106.   IK = 100L*(N-49L) + IT + L;
  107.   *mon = (int)JT; *day = (int)K; *yea = (int)IK;
  108. #else
  109.   double tim;
  110.  
  111.   revjul(JD, JD >= 2299171.0 /* October 15, 1582 */, mon, day, yea, &tim);
  112. #endif
  113. }
  114.  
  115.  
  116. /* This is a subprocedure of CastChart(). Once we have the chart parameters, */
  117. /* calculate a few important things related to the date, i.e. the Greenwich  */
  118. /* time, the Julian day and fractional part of the day, the offset to the    */
  119. /* sidereal, and a couple of other things.                                   */
  120.  
  121. real ProcessInput(fDate)
  122. bool fDate;
  123. {
  124.   real Off, Ln;
  125.  
  126.   TT = RSgn(TT)*RFloor(RAbs(TT))+RFract(RAbs(TT))*100.0/60.0 +
  127.     (DecToDeg(ZZ) - DecToDeg(SS));
  128.   OO = DecToDeg(OO);
  129.   AA = Min(AA, 89.9999);        /* Make sure the chart isn't being cast */
  130.   AA = Max(AA, -89.9999);       /* on the precise north or south pole.  */
  131.   AA = RFromD(DecToDeg(AA));
  132.  
  133.   /* if parameter 'fDate' isn't set, then we can assume that the true time */
  134.   /* has already been determined (as in a -rm switch time midpoint chart). */
  135.  
  136.   if (fDate) {
  137.     is.JD = (real)MdyToJulian(MM, DD, YY);
  138.     if (!us.fProgress || us.fSolarArc)
  139.       T = (is.JD + TT/24.0 - 2415020.5) / 36525.0;
  140.     else {
  141.       /* Determine actual time that a progressed chart is to be cast for. */
  142.  
  143.       T = ((is.JD + TT/24.0 + (is.JDp - (is.JD + TT/24.0)) / us.rProgDay) -
  144.         2415020.5) / 36525.0;
  145.     }
  146.   }
  147.  
  148.   /* Compute angle that the ecliptic is inclined to the Celestial Equator */
  149.   OB = RFromD(23.452294-0.0130125*T);
  150.  
  151.   Ln = Mod((933060-6962911*T+7.5*T*T)/3600.0);    /* Mean lunar node */
  152.   Off = (259205536.0*T+2013816.0)/3600.0;         /* Mean Sun        */
  153.   Off = 17.23*RSin(RFromD(Ln))+1.27*RSin(RFromD(Off))-(5025.64+1.11*T)*T;
  154.   Off = (Off-84038.27)/3600.0;
  155.   is.rSid = (us.fSiderial ? Off : 0.0) + us.rZodiacOffset;
  156.   return Off;
  157. }
  158.  
  159.  
  160. /* Convert polar to rectangular coordinates. */
  161.  
  162. void PolToRec(A, R, X, Y)
  163. real A, R, *X, *Y;
  164. {
  165.   if (A == 0.0)
  166.     A = rSmall;
  167.   *X = R*RCos(A);
  168.   *Y = R*RSin(A);
  169. }
  170.  
  171.  
  172. /* Convert rectangular to polar coordinates. */
  173.  
  174. void RecToPol(X, Y, A, R)
  175. real X, Y, *A, *R;
  176. {
  177.   if (Y == 0.0)
  178.     Y = rSmall;
  179.   *R = RSqr(X*X + Y*Y);
  180.   *A = Angle(X, Y);
  181. }
  182.  
  183.  
  184. /* Convert rectangular to spherical coordinates. */
  185.  
  186. real RecToSph(B, L, O)
  187. real B, L, O;
  188. {
  189.   real R, Q, G, X, Y, A;
  190.  
  191.   A = B; R = 1.0;
  192.   PolToRec(A, R, &X, &Y);
  193.   Q = Y; R = X; A = L;
  194.   PolToRec(A, R, &X, &Y);
  195.   G = X; X = Y; Y = Q;
  196.   RecToPol(X, Y, &A, &R);
  197.   A += O;
  198.   PolToRec(A, R, &X, &Y);
  199.   Q = RAsin(Y);
  200.   Y = X; X = G;
  201.   RecToPol(X, Y, &A, &R);
  202.   if (A < 0.0)
  203.     A += 2*rPi;
  204.   G = A;
  205.   return G;  /* We only ever care about and return one of the coordinates. */
  206. }
  207.  
  208.  
  209. /* Do a coordinate transformation: Given a longitude and latitude value,    */
  210. /* return the new longitude and latitude values that the same location      */
  211. /* would have, were the equator tilted by a specified number of degrees.    */
  212. /* In other words, do a pole shift! This is used to convert among ecliptic, */
  213. /* equatorial, and local coordinates, each of which have zero declination   */
  214. /* in different planes. In other words, take into account the Earth's axis. */
  215.  
  216. void CoorXform(azi, alt, tilt)
  217. real *azi, *alt, tilt;
  218. {
  219.   real x, y, a1, l1;
  220.   real sinalt, cosalt, sinazi, sintilt, costilt;
  221.  
  222.   sinalt = RSin(*alt); cosalt = RCos(*alt); sinazi = RSin(*azi);
  223.   sintilt = RSin(tilt); costilt = RCos(tilt);
  224.  
  225.   x = cosalt * sinazi * costilt;
  226.   y = sinalt * sintilt;
  227.   x -= y;
  228.   a1 = cosalt;
  229.   y = cosalt * RCos(*azi);
  230.   l1 = Angle(y, x);
  231.   a1 = a1 * sinazi * sintilt + sinalt * costilt;
  232.   a1 = RAsin(a1);
  233.   *azi = l1; *alt = a1;
  234. }
  235.  
  236.  
  237. /* This is another subprocedure of CastChart(). Calculate a few variables */
  238. /* corresponding to the chart parameters that are used later on. The      */
  239. /* astrological vertex (object number twenty) is also calculated here.    */
  240.  
  241. void ComputeVariables(vtx)
  242. real *vtx;
  243. {
  244.   real R, R2, B, L, O, G, X, Y, A;
  245.  
  246.   RA = RFromD(Mod((6.6460656+2400.0513*T+2.58E-5*T*T+TT)*15.0-OO));
  247.   R2 = RA; O = -OB; B = AA; A = R2; R = 1.0;
  248.   PolToRec(A, R, &X, &Y);
  249.   X *= RCos(O);
  250.   RecToPol(X, Y, &A, &R);
  251.   MC = Mod(is.rSid + DFromR(A));              /* Midheaven */
  252. #if FALSE
  253.   L = R2;
  254.   G = RecToSph(B, L, O);
  255.   Asc = Mod(is.rSid + Mod(G+rPiHalf));        /* Ascendant */
  256. #endif
  257.   L= R2+rPi; B = rPiHalf-RAbs(B);
  258.   if (AA < 0.0)
  259.     B = -B;
  260.   G = RecToSph(B, L, O);
  261.   *vtx = Mod(is.rSid + DFromR(G+rPiHalf));    /* Vertex */
  262. }
  263.  
  264.  
  265. /*
  266. ******************************************************************************
  267. ** House Cusp Calculations.
  268. ******************************************************************************
  269. */
  270.  
  271. /* The following three functions calculate the Midheaven, Ascendant, and  */
  272. /* East Point of the chart in question, based on time and location. The   */
  273. /* first two are also used in some of the house cusp calculation routines */
  274. /* as a quick way to get the 10th and 1st cusps. The East Point object is */
  275. /* technically defined as the Ascendant's position at zero latitude.      */
  276.  
  277. real CuspMidheaven()
  278. {
  279.   real MC;
  280.  
  281.   MC = RAtn(RTan(RA)/RCos(OB));
  282.   if (MC < 0.0)
  283.     MC += rPi;
  284.   if (RA > rPi)
  285.     MC += rPi;
  286.   return Mod(DFromR(MC)+is.rSid);
  287. }
  288.  
  289. real CuspAscendant()
  290. {
  291.   real Asc;
  292.  
  293.   Asc = Angle(-RSin(RA)*RCos(OB)-RTan(AA)*RSin(OB), RCos(RA));
  294.   return Mod(DFromR(Asc)+is.rSid);
  295. }
  296.  
  297. real CuspEastPoint()
  298. {
  299.   real EP;
  300.  
  301.   EP = Angle(-RSin(RA)*RCos(OB), RCos(RA));
  302.   return Mod(DFromR(EP)+is.rSid);
  303. }
  304.  
  305.  
  306. /* These are various different algorithms for calculating the house cusps: */
  307.  
  308. real CuspPlacidus(deg, FF, fNeg)
  309. real deg, FF;
  310. bool fNeg;
  311. {
  312.   real LO, R1, XS, X;
  313.   int i;
  314.  
  315.   R1 = RA+RFromD(deg);
  316.   X = fNeg ? 1.0 : -1.0;
  317.   /* Looping 10 times is arbitrary, but it's what other programs do. */
  318.   for (i = 1; i <= 10; i++) {
  319.  
  320.     /* This formula works except at 0 latitude (AA == 0.0). */
  321.  
  322.     XS = X*RSin(R1)*RTan(OB)*RTan(AA == 0.0 ? 0.0001 : AA);
  323.     XS = RAcos(XS);
  324.     if (XS < 0.0)
  325.       XS += rPi;
  326.     R1 = RA + (fNeg ? rPi-(XS/FF) : (XS/FF));
  327.   }
  328.   LO = RAtn(RTan(R1)/RCos(OB));
  329.   if (LO < 0.0)
  330.     LO += rPi;
  331.   if (RSin(R1) < 0.0)
  332.     LO += rPi;
  333.   return DFromR(LO);
  334. }
  335.  
  336. void HousePlacidus()
  337. {
  338.   int i;
  339.  
  340.   house[1] = Mod(Asc-is.rSid);
  341.   house[4] = Mod(MC+rDegHalf-is.rSid);
  342.   house[5] = CuspPlacidus(30.0, 3.0, fFalse) + rDegHalf;
  343.   house[6] = CuspPlacidus(60.0, 1.5, fFalse) + rDegHalf;
  344.   house[2] = CuspPlacidus(120.0, 1.5, fTrue);
  345.   house[3] = CuspPlacidus(150.0, 3.0, fTrue);
  346.   for (i = 1; i <= cSign; i++) {
  347.     if (i <= 6)
  348.       house[i] = Mod(house[i]+is.rSid);
  349.     else
  350.       house[i] = Mod(house[i-6]+rDegHalf);
  351.   }
  352. }
  353.  
  354. void HouseKoch()
  355. {
  356.   real A1, A2, A3, KN, D, X;
  357.   int i;
  358.  
  359.   A1 = RSin(RA)*RTan(AA)*RTan(OB);
  360.   A1 = RAsin(A1);
  361.   for (i = 1; i <= cSign; i++) {
  362.     D = Mod(60.0+30.0*(real)i);
  363.     A2 = D/rDegQuad-1.0; KN = 1.0;
  364.     if (D >= rDegHalf) {
  365.       KN = -1.0;
  366.       A2 = D/rDegQuad-3.0;
  367.     }
  368.     A3 = RFromD(Mod(DFromR(RA)+D+A2*DFromR(A1)));
  369.     X = Angle(RCos(A3)*RCos(OB)-KN*RTan(AA)*RSin(OB), RSin(A3));
  370.     house[i] = Mod(DFromR(X)+is.rSid);
  371.   }
  372. }
  373.  
  374. void HouseEqual()
  375. {
  376.   int i;
  377.  
  378.   for (i = 1; i <= cSign; i++)
  379.     house[i] = Mod(Asc-30.0+30.0*(real)i);
  380. }
  381.  
  382. void HouseCampanus()
  383. {
  384.   real KO, DN, X;
  385.   int i;
  386.  
  387.   for (i = 1; i <= cSign; i++) {
  388.     KO = RFromD(60.000001+30.0*(real)i);
  389.     DN = RAtn(RTan(KO)*RCos(AA));
  390.     if (DN < 0.0)
  391.       DN += rPi;
  392.     if (RSin(KO) < 0.0)
  393.       DN += rPi;
  394.     X = Angle(RCos(RA+DN)*RCos(OB)-RSin(DN)*RTan(AA)*RSin(OB), RSin(RA+DN));
  395.     house[i] = Mod(DFromR(X)+is.rSid);
  396.   }
  397. }
  398.  
  399. void HouseMeridian()
  400. {
  401.   real D, X;
  402.   int i;
  403.  
  404.   for (i = 1; i <= cSign; i++) {
  405.     D = RFromD(60.0+30.0*(real)i);
  406.     X = Angle(RCos(RA+D)*RCos(OB), RSin(RA+D));
  407.     house[i] = Mod(DFromR(X)+is.rSid);
  408.   }
  409. }
  410.  
  411. void HouseRegiomontanus()
  412. {
  413.   real D, X;
  414.   int i;
  415.  
  416.   for (i = 1; i <= cSign; i++) {
  417.     D = RFromD(60.0+30.0*(real)i);
  418.     X = Angle(RCos(RA+D)*RCos(OB)-RSin(D)*RTan(AA)*RSin(OB), RSin(RA+D));
  419.     house[i] = Mod(DFromR(X)+is.rSid);
  420.   }
  421. }
  422.  
  423. void HousePorphyry()
  424. {
  425.   real X, Y;
  426.   int i;
  427.  
  428.   X = Asc-MC;
  429.   if (X < 0.0)
  430.     X += rDegMax;
  431.   Y = X/3.0;
  432.   for (i = 1; i <= 2; i++)
  433.     house[i+4] = Mod(rDegHalf+MC+i*Y);
  434.   X = Mod(rDegHalf+MC)-Asc;
  435.   if (X < 0.0)
  436.     X += rDegMax;
  437.   house[1]=Asc;
  438.   Y = X/3.0;
  439.   for (i = 1; i <= 3; i++)
  440.     house[i+1] = Mod(Asc+i*Y);
  441.   for (i = 1; i <= 6; i++)
  442.     house[i+6] = Mod(house[i]+rDegHalf);
  443. }
  444.  
  445. void HouseMorinus()
  446. {
  447.   real D, X;
  448.   int i;
  449.  
  450.   for (i = 1; i <= cSign; i++) {
  451.     D = RFromD(60.0+30.0*(real)i);
  452.     X = Angle(RCos(RA+D), RSin(RA+D)*RCos(OB));
  453.     house[i] = Mod(DFromR(X)+is.rSid);
  454.   }
  455. }
  456.  
  457. real CuspTopocentric(deg)
  458. real deg;
  459. {
  460.   real OA, X, LO;
  461.  
  462.   OA = ModRad(RA+RFromD(deg));
  463.   X = RAtn(RTan(AA)/RCos(OA));
  464.   LO = RAtn(RCos(X)*RTan(OA)/RCos(X+OB));
  465.   if (LO < 0.0)
  466.     LO += rPi;
  467.   if (RSin(OA) < 0.0)
  468.     LO += rPi;
  469.   return LO;
  470. }
  471.  
  472. void HouseTopocentric()
  473. {
  474.   real TL, P1, P2, LT;
  475.   int i;
  476.  
  477.   house[4] = ModRad(RFromD(MC+rDegHalf-is.rSid));
  478.   TL = RTan(AA); P1 = RAtn(TL/3.0); P2 = RAtn(TL/1.5); LT = AA;
  479.   AA = P1; house[5] = CuspTopocentric(30.0) + rPi;
  480.   AA = P2; house[6] = CuspTopocentric(60.0) + rPi;
  481.   AA = LT; house[1] = CuspTopocentric(90.0);
  482.   AA = P2; house[2] = CuspTopocentric(120.0);
  483.   AA = P1; house[3] = CuspTopocentric(150.0);
  484.   AA = LT;
  485.   for (i = 1; i <= 6; i++) {
  486.     house[i] = Mod(DFromR(house[i])+is.rSid);
  487.     house[i+6] = Mod(house[i]+rDegHalf);
  488.   }
  489. }
  490.  
  491.  
  492. /*
  493. ******************************************************************************
  494. ** Planetary Position Calculations.
  495. ******************************************************************************
  496. */
  497.  
  498. /* Given three values, return them combined as the coefficients of a */
  499. /* quadratic equation as a function of the chart time.               */
  500.  
  501. real ReadThree(r0, r1, r2)
  502. real r0, r1, r2;
  503. {
  504.   return RFromD(r0 + r1*T + r2*T*T);
  505. }
  506.  
  507.  
  508. /* Another coordinate transformation. This is used by the ComputePlanets() */
  509. /* procedure to rotate rectangular coordinates by a certain amount.        */
  510.  
  511. void RecToSph2(AP, AN, IN, X, Y, G)
  512. real AP, AN, IN, *X, *Y, *G;
  513. {
  514.   real R, D, A;
  515.  
  516.   RecToPol(*X, *Y, &A, &R); A += AP; PolToRec(A, R, X, Y);
  517.   D = *X; *X = *Y; *Y = 0.0; RecToPol(*X, *Y, &A, &R);
  518.   A += IN; PolToRec(A, R, X, Y);
  519.   *G = *Y; *Y = *X; *X = D; RecToPol(*X, *Y, &A, &R); A += AN;
  520.   if (A < 0.0)
  521.     A += 2.0*rPi;
  522.   PolToRec(A, R, X, Y);
  523. }
  524.  
  525.  
  526. /* Calculate some harmonic delta error correction factors to add onto the */
  527. /* coordinates of Jupiter through Pluto, for better accuracy.             */
  528.  
  529. void ErrorCorrect(ind, x, y, z)
  530. int ind;
  531. real *x, *y, *z;
  532. {
  533.   real U, V, W, A, S0, T0[4], FAR *pr;
  534.   int IK, IJ, irError;
  535.  
  536.   irError = errorcount[ind-oJup];
  537.   pr = (lpreal)&errordata[erroroffset[ind-oJup]];
  538.   for (IK = 1; IK <= 3; IK++) {
  539.     if (ind == oJup && IK == 3) {
  540.       T0[3] = 0.0;
  541.       break;
  542.     }
  543.     if (IK == 3)
  544.       irError--;
  545.     S0 = ReadThree(pr[0], pr[1], pr[2]); pr += 3;
  546.     A = 0.0;
  547.     for (IJ = 1; IJ <= irError; IJ++) {
  548.       U = *pr++; V = *pr++; W = *pr++;
  549.       A += RFromD(U)*RCos((V*T+W)*rPi/rDegHalf);
  550.     }
  551.     T0[IK] = DFromR(S0+A);
  552.   }
  553.   *x += T0[2]; *y += T0[1]; *z += T0[3];
  554. }
  555.  
  556.  
  557. /* Another subprocedure of the ComputePlanets() routine. Convert the final */
  558. /* rectangular coordinates of a planet to zodiac position and declination. */
  559.  
  560. void ProcessPlanet(ind, aber)
  561. int ind;
  562. real aber;
  563. {
  564.   real ang, rad;
  565.  
  566.   RecToPol(spacex[ind], spacey[ind], &ang, &rad);
  567.   planet[ind] = Mod(DFromR(ang) /*+ NU*/ - aber + is.rSid);
  568.   RecToPol(rad, spacez[ind], &ang, &rad);
  569.   if (us.objCenter == oSun && ind == oSun)
  570.     ang = 0.0;
  571.   ang = DFromR(ang);
  572.   while (ang > rDegQuad)    /* Ensure declination is from -90..+90 degrees. */
  573.     ang -= rDegHalf;
  574.   while (ang < -rDegQuad)
  575.     ang += rDegHalf;
  576.   planetalt[ind] = ang;
  577. }
  578.  
  579.  
  580. /* This is probably the heart of the whole program of Astrolog. Calculate  */
  581. /* the position of each body that orbits the Sun. A heliocentric chart is  */
  582. /* most natural; extra calculation is needed to have other central bodies. */
  583.  
  584. void ComputePlanets()
  585. {
  586.   real helio[oNorm+1], helioalt[oNorm+1], helioret[oNorm+1],
  587.     heliox[oNorm+1], helioy[oNorm+1], helioz[oNorm+1];
  588.   real aber = 0.0, AU, E, EA, E1, M, XW, YW, AP, AN, IN, X, Y, G, XS, YS, ZS;
  589.   int ind = oSun, i;
  590.   OE *poe;
  591.  
  592.   while (ind <= (us.fUranian ? oNorm : cPlanet)) {
  593.     if (ignore[ind] && ind > oSun)
  594.       goto LNextPlanet;
  595.     poe = &rgoe[IoeFromObj(ind)];
  596.  
  597.     EA = M = ModRad(ReadThree(poe->ma0, poe->ma1, poe->ma2));
  598.     E = DFromR(ReadThree(poe->ec0, poe->ec1, poe->ec2));
  599.     for (i = 1; i <= 5; i++)
  600.       EA = M+E*RSin(EA);            /* Solve Kepler's equation */
  601.     AU = poe->sma;                  /* Semi-major axis         */
  602.     E1 = 0.01720209/(pow(AU, 1.5)*
  603.       (1.0-E*RCos(EA)));                     /* Begin velocity coordinates */
  604.     XW = -AU*E1*RSin(EA);                    /* Perifocal coordinates      */
  605.     YW = AU*E1*pow(1.0-E*E,0.5)*RCos(EA);
  606.     AP = ReadThree(poe->ap0, poe->ap1, poe->ap2);
  607.     AN = ReadThree(poe->an0, poe->an1, poe->an2);
  608.     IN = ReadThree(poe->in0, poe->in1, poe->in2); /* Calculate inclination  */
  609.     X = XW; Y = YW;
  610.     RecToSph2(AP, AN, IN, &X, &Y, &G);  /* Rotate velocity coords */
  611.     heliox[ind] = X; helioy[ind] = Y;
  612.     helioz[ind] = G;                    /* Helio ecliptic rectangtular */
  613.     X = AU*(RCos(EA)-E);                /* Perifocal coordinates for        */
  614.     Y = AU*RSin(EA)*pow(1.0-E*E,0.5);   /* rectangular position coordinates */
  615.     RecToSph2(AP, AN, IN, &X, &Y, &G);  /* Rotate for rectangular */
  616.     XS = X; YS = Y; ZS = G;             /* position coordinates   */
  617.     if (FBetween(ind, oJup, oPlu))
  618.       ErrorCorrect(ind, &XS, &YS, &ZS);
  619.     ret[ind] =                                        /* Helio daily motion */
  620.       (XS*helioy[ind]-YS*heliox[ind])/(XS*XS+YS*YS);
  621.     spacex[ind] = XS; spacey[ind] = YS; spacez[ind] = ZS;
  622.     ProcessPlanet(ind, 0.0);
  623. LNextPlanet:
  624.     ind += (ind == oSun ? 2 : (ind != cPlanet ? 1 : uranLo-cPlanet));
  625.   }
  626.   spacex[0] = spacey[0] = spacez[0] = 0.0;
  627.  
  628.   if (!us.objCenter) {
  629.     if (us.fVelocity)
  630.       for (i = 0; i <= oNorm; i++)    /* Use relative velocity */
  631.         ret[i] = RFromD(1.0);         /* if -v0 is in effect.  */
  632.     return;
  633.   }
  634. #endif /* MATRIX */
  635.  
  636.   /* A second loop is needed for geocentric charts or central bodies other */
  637.   /* than the Sun. For example, we can't find the position of Mercury in   */
  638.   /* relation to Pluto until we know the position of Pluto in relation to  */
  639.   /* the Sun, and since Mercury is calculated first, another pass needed.  */
  640.  
  641.   ind = us.objCenter;
  642.   for (i = 0; i <= oNorm; i++) {
  643.     helio[i]    = planet[i];
  644.     helioalt[i] = planetalt[i];
  645.     helioret[i] = ret[i];
  646.     if (i != oMoo && i != ind) {
  647.       spacex[i] -= spacex[ind];
  648.       spacey[i] -= spacey[ind];
  649.       spacez[i] -= spacez[ind];
  650.     }
  651.   }
  652.   spacex[ind] = spacey[ind] = spacez[ind] = 0.0;
  653.   SwapR(&spacex[0], &spacex[ind]);
  654.   SwapR(&spacey[0], &spacey[ind]);    /* Do some swapping - we want   */
  655.   SwapR(&spacez[0], &spacez[ind]);    /* the central body to be in    */
  656.   SwapR(&spacex[1], &spacex[ind]);    /* object position number zero. */
  657.   SwapR(&spacey[1], &spacey[ind]);
  658.   SwapR(&spacez[1], &spacez[ind]);
  659.   for (i = 1; i <= (us.fUranian ? oNorm : cPlanet);
  660.       i += (i == 1 ? 2 : (i != cPlanet ? 1 : uranLo-cPlanet))) {
  661.     if (ignore[i] && i > oSun)
  662.       continue;
  663.     XS = spacex[i]; YS = spacey[i]; ZS = spacez[i];
  664.     if (ind != oSun || i != oSun)
  665.       ret[i] = (XS*(helioy[i]-helioy[ind])-YS*(heliox[i]-heliox[ind]))/
  666.         (XS*XS+YS*YS);
  667.     if (ind == oSun)
  668.       aber = 0.0057756*RSqr(XS*XS+YS*YS+ZS*ZS)*DFromR(ret[i]); /* Aberration */
  669.     ProcessPlanet(i, aber);
  670.     if (us.fVelocity)                         /* Use relative velocity */
  671.       ret[i] = RFromD(ret[i]/helioret[i]);    /* if -v0 is in effect   */
  672.   }
  673. }
  674.  
  675.  
  676. #ifdef MATRIX
  677. /*
  678. ******************************************************************************
  679. ** Lunar Position Calculations
  680. ******************************************************************************
  681. */
  682.  
  683. /* Calculate the position and declination of the Moon, and the Moon's North  */
  684. /* Node. This has to be done separately from the other planets, because they */
  685. /* all orbit the Sun, while the Moon orbits the Earth.                       */
  686.  
  687. void ComputeLunar(moonlo, moonla, nodelo, nodela)
  688. real *moonlo, *moonla, *nodelo, *nodela;
  689. {
  690.   real LL, G, N, G1, D, L, ML, L1, MB, T1, Y, M = 3600.0;
  691.  
  692.   LL = 973563.0+1732564379.0*T-4.0*T*T;  /* Compute mean lunar longitude     */
  693.   G = 1012395.0+6189.0*T;                /* Sun's mean longitude of perigee  */
  694.   N = 933060.0-6962911.0*T+7.5*T*T;      /* Compute mean lunar node          */
  695.   G1 = 1203586.0+14648523.0*T-37.0*T*T;  /* Mean longitude of lunar perigee  */
  696.   D = 1262655.0+1602961611.0*T-5.0*T*T;  /* Mean elongation of Moon from Sun */
  697.   L = (LL-G1)/M; L1 = ((LL-D)-G)/M;      /* Some auxiliary angles            */
  698.   T1 = (LL-N)/M; D = D/M; Y = 2.0*D;
  699.  
  700.   /* Compute Moon's perturbations. */
  701.  
  702.   ML = 22639.6*RSinD(L) - 4586.4*RSinD(L-Y) + 2369.9*RSinD(Y) +
  703.     769.0*RSinD(2.0*L) - 669.0*RSinD(L1) - 411.6*RSinD(2.0*T1) -
  704.     212.0*RSinD(2.0*L-Y) - 206.0*RSinD(L+L1-Y);
  705.   ML += 192.0*RSinD(L+Y) - 165.0*RSinD(L1-Y) + 148.0*RSinD(L-L1) -
  706.     125.0*RSinD(D) - 110.0*RSinD(L+L1) - 55.0*RSinD(2.0*T1-Y) -
  707.     45.0*RSinD(L+2.0*T1) + 40.0*RSinD(L-2.0*T1);
  708.  
  709.   *moonlo = G = Mod((LL+ML)/M+is.rSid);    /* Lunar longitude */
  710.  
  711.   /* Compute lunar latitude. */
  712.  
  713.   MB = 18461.5*RSinD(T1) + 1010.0*RSinD(L+T1) - 999.0*RSinD(T1-L) -
  714.     624.0*RSinD(T1-Y) + 199.0*RSinD(T1+Y-L) - 167.0*RSinD(L+T1-Y);
  715.   MB += 117.0*RSinD(T1+Y) + 62.0*RSinD(2.0*L+T1) -
  716.     33.0*RSinD(T1-Y-L) - 32.0*RSinD(T1-2.0*L) - 30.0*RSinD(L1+T1-Y);
  717.   *moonla = MB =
  718.     RSgn(MB)*((RAbs(MB)/M)/rDegMax-RFloor((RAbs(MB)/M)/rDegMax))*rDegMax;
  719.  
  720.   /* Compute position of the North Lunar Node, either True or Mean. */
  721.  
  722.   if (us.fTrueNode)
  723.     N = N+5392.0*RSinD(2.0*T1-Y)-541.0*RSinD(L1)-442.0*RSinD(Y)+
  724.       423.0*RSinD(2.0*T1)-291.0*RSinD(2.0*L-2.0*T1);
  725.   *nodelo = Mod(N/M+is.rSid);
  726.   *nodela = 0.0;
  727. }
  728. #endif /* MATRIX */
  729.  
  730. /* matrix.c */
  731.